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We investigate the properties of hard core Bosons in harmonic traps over a wide range of densities. 
Bose-Einstein condensation is formulated using the one-body Density Matrix (OBDM) which is 
equally valid at low and high densities. The OBDM is calculated using diffusion Monte Carlo 
methods and it is diagonalized to obtain the "natural" single particle orbitals and their occupation, 
including the condensate fraction. At low Boson density, na 3 < 10~ 5 , where n = N/V and a is 
the hard core diameter, the condensate is localized at the center of the trap. As na 3 increases, the 
condensate moves to the edges of the trap. At high density it is localized at the edges of the trap. 
At na 3 < 10~ 4 the Gross-Pitaevskii theory of the condensate describes the whole system within 
1%. At na 3 « 10~ 3 corrections are 3% to the GP energy but 30% to the Bogoliubov prediction of 
the condensate depletion. At na 3 > 10 -2 , mean field theory fails. At na 3 > 0.1, the Bosons behave 
more like a liquid 4 He droplet than a trapped Boson gas. 



I. INTRODUCTION 

Bose-Einstein condensation (BEC) has been a topic of 
fundamental interest since it was first predicted by Ein- 
stein in 1924 [| . He showed that as a consequence of Bose 
statistics Q a macroscopic fraction, Nq/N, of the atoms 
in an ideal Bosegas can condense into a single quantum 
state. London 0, Q] postulated that superfluidity in liq- 
uid 4 He was a consequence of a transition to BEC. But 
liquid 4 He is a strongly interacting, dense Bose liquid and 
the connection between BEC in an ideal gas and super- 
fluidity was not at all clear |{| . Similarly, the many-body 
correlation effects induced by the inter-boson interaction 
significantly reduce the condensate fraction even at zero 
temperature [(J Q • Modern direct measurements of 
BEC in liquid 4 He find only 7.25% of the liquid in the 
condensate at T = OK. 

The theoretical framework for treating an interacting 
Bose gas was initiated in 1947 by Bogoliubov 0- He 
developed a perturbation expansion valid for low den- 
sity and weak interaction, na 3 -C 1 (where n is the 
number density N/V and a is the hard core diameter 
of the Bosons), and small depletion of the condensate, 
(N — Nn)/N <c 1. About a decade later, Onsager and 
Penrose and Lowdin [Tl| formulated a definition of 
BEC in terms of the eigenvalues and eigenvectors (nat- 
ural orbitals) of the one-body density matrix (OBDM). 
An orbital with macroscopic occupation arising from di- 
agonalization of the OBDM is defined as the "condensate 
wave- function" or order parameter. This formulation al- 
lows direct access to condensate properties at arbitrary 
density and does not require a large condensate fraction. 
The work in this paper is based on the OBDM formu- 
lation of BEC which is rigorously valid for a strongly 
interacting system |5j. 

In 1995, experiments in weakly interacting dilute va- 
pors of the alkali atoms 87 Rb, 23 Na and 7 Li in magnetic 
traps provided direct evidence of a clear transition from 
a thermally distributed cloud to macroscopic occupation 
of a single quantum state 0, 0, . This long awaited 
direct realization of BEC spawned a dramatic renewal 



of interest in Bose systems and BEC. Since the densities 
in these experiments were low (typical number densities 
were 10 12 cm~ 3 and na 3 ~ 10~ 6 where a is the s-wave 
scattering length of the atoms), almost all of the theo- 
retical activity has focused on the weakly interacting gas 
limit and the Gross-Pitaevskii (GP) equation [Tfij |. The 
GP equation provides an excellent mean field description 
of the condensate at low density. This is a valid descrip- 
tion of the whole Bose gas in the dilute limit, na 3 <C 1, 
where most of the atoms are in the condensate. How- 
ever, it is inaccurate for strongly interacting systems in 
which the condensate fraction is significantly depleted 
by quantum fluctuations. Since the experiments in 1995, 
only a handful of studies have attempted to consider the 
properties of BEC beyond the dilute regime and the GP 
description of the condensate PbSl ITtI flip. \l$L I20L l2~p. I22I 

EiMmmmmmiM mr iimim sua m . 

Most of this relatively small body of work rely on mod- 
ified forms of the GP equation which incorporate higher 
terms in the Bogoliubov expansion that include effects 
of atoms outside the condensate within a local density 
approximation. Unfortunately, the condensate fraction 
and distribution in the trap calculated by such methods 
become inaccurate as the density becomes greater than 
na 3 > 1CT 3 HI. 

It has recently become possible to study Bose systems 
with tunable interactions [H HHU IS El El for which 
densities of up to na 3 ss 1 are obtainable. Specifically, 
85 Rb at densities in the range na 3 « 1CP 3 — 10 _1 has been 
investigated. BEC in metastable helium isotopes ELIiol 
l4f| | with na 3 ~ 10 -4 and in atomic hydrogen |47| with 
na 3 10~ 5 , are also higher density Bose gases. This 
makes the study of BEC and the role of interactions in 
trapped Bose gases over a wide range of densities of direct 
interest to experiment. 

The chief purpose of this work is to go beyond the 
dilute limit, to test the limits of the GP equation and re- 
lated mean field approximations and to explore the zero 
temperature properties of trapped hard core Bosons as 
na 3 increases from the dilute limit to the dense regime 
corresponding to liquid 4 He, and beyond. The range of 
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FIG. 1: Range of system densities considered in this work 
expressed in terms of na s = Na 3 /V, the ratio of the volume 
occupied by N hard core particles with diameter a to the total 
volume of the system, V. 



densities investigated here is displayed in Fig. ^ We in- 
crease the density by increasing both N and the hard core 
diameter, a, up to the value na 3 ~ 0.21 which describes 
liquid 4 He at SVP when the 4 He atoms are represented 
by hard spheres of diameter a = 2.203A |^. The ground 
state energy, E, the total density distribution, and the 
OBDM are evaluated using diffusion Monte Carlo (DMC) 
methods. 

Specifically, we compare the ground state energy of 
the whole trapped gas calculated using DMC, Edmc-, 
with the usual energy of the condensate calculated using 
the GP equation, Eqp- As density increases, Edmc and 
Eqp begin to differ. For example, at na 3 — 10~ 3 , we find 
{Edmc - E GP ) / E GP = 3%. Modified GP equations pro- 
vide a mean field description of the atoms above the con- 
densate. The dependence of Eumc~^GP on the number 
of trapped Bosons, N , and on the scattering length, a, 
follows the predictions of the Modified GP equation re- 
markably well up to high densities, na 3 » 5 x 10" 2 . This 
suggests that the difference Edmc — Evmc can be at- 
tributed to the atoms above the condensate. However, 
the energy is not as sensitive to approximations as some 
other properties. 

We compare the condensate fraction obtained using 
the rigorous OBDM-DMC method with predictions of 
the Bogoliubov theory. The two agree within 1% for 
na 3 < 10~ 4 . At higher densities, the Bogoliubov theory 
significantly underestimates the depletion of the conden- 
sate, by 25% at na 3 w 2 x 10~ 2 . We evaluate the con- 
densate density distribution in the trap. At low density, 
the condensate is localized at the center of the trap as 
usually found ^5|. At higher density (na 3 10~ 2 ), the 
condensate is spaced over several trap lengths and the 
condensate and total density have similar distributions. 
Also, at higher densities (na 3 >2x 10~ 2 ), oscillations in 
the total density distribution appear which are not found 
in mean field theories. There are no corresponding oscil- 
lations in the condensate density distribution. At high 
density (na 3 > 0.10), the condensate is localized at the 
edges of the trap (large r/a) where the total boson den- 
sity is low. At high density, the trapped Bosons resemble 
liquid 4 He droplets HI Hfl HU . 

We also compare the present DMC results with our 
earlier variational Monte Carlo (VMC) values [5(J. We 



find that the VMC and DMC energies agree well at all 
densities and that the ground state energy is not very 
sensitive to the trial variational wave-function. However, 
the OBDM and the condensate fraction is very sensitive 
to the trial wave function at higher densities. An accurate 
initial trial function is needed to get reliable condensate 
fractions even in the DMC formulation. 

Monte Carlo methods are usually applied to dense sys- 
tems such as liquid and solid 4 He [a, Q . Recently Giorgini 
et al. have evaluated the energy and condensate frac- 
tion of the uniform Bose gas over a wide density range, 
10~ 6 < na 3 < 10" 3 0. Griiter et al. HJ, have evalu- 
ated the critical temperature, T c , for BEC in a Bose gas 
using path integral Monte Carlo (PIMC) methods. They 
find T c is increased above the ideal Bose gas value by 
interaction in the dilute range. This increase is observed 
in dilute concentrations of 4 He in Vycor At liquid 
4 He densities, T c is decreased by interaction |53l IH^j. 

Krauth [16] first applied QMC to BEC in a trap us- 
ing PIMC methods. For 10, 000 hard sphere Bosons in a 
spherical trap with a ratio of hard core diameter to trap 
length, a/a ho = 4.3 x 10 -3 (na 3 w 10~ 4 ), he found that 
condensate was concentrated at the center of the trap 
while the uncondensed atoms were spread over a wide 
range and well described by a classical Bose gas. Holz- 
mann et al. j2^| made a direct comparison of PIMC and 
Hartree-Fock calculations for a dilute gas of hard spheres 
in a trap with a/a% — 4.3 x 10 -3 . For temperatures near 
T c , they found Nq was greater in PIMC. The increase in 
Nq with exact representation of the interaction effects 
is consistent with the corresponding increase in T c with 
interaction in the uniform Bose gas. 

Recently, QMC methods have been successfully ap- 
plied to the study of highly inhomogeneous Bose systems. 
Astrakharchik et al. used DMC to study BEC and super- 
fluidity in a Bose gas with disorder at zero temperature 
[3ll ] . They find an intriguing decoupling of the superfluid 
and condensate fractions for strong disorder. Studies of 
superfluid 4 He with a free surface [27], l33| found the local 
condensate fraction peaks (no ~ 0.95% |33^~) in the dilute 
region just inside the liquid- vacuum interface. Blume j3|| 
and Astrakharchik and Giorgini 37] have examined the 
transition from the three-dimensional to the quasi one- 
dimensional regime for Bosons in highly elongated cigar- 
shaped traps. They confirm that the Bose gas undergoes 
"fermionization" in the quasi 1-D regime. 

In section [H] we describe the theoretical framework 
and computational methods used. Section IIIII contains 
the present results. In section IIVI the chief results are 
reviewed and discussed. 



II. METHODS 

We consider N Bosons of mass m confined in an ex- 
ternal trapping potential, V ex t('r), and interacting via a 
two-body potential Vint ( r i , r 2) • The Hamiltonian for this 
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system is: 

N , ,2 



N , 2 v AT 



Here, 



V ext (r) = -mu;l r 2 , 



(2) 



where u> 2 is the characteristic trap frequency. Interac- 
tions are modelled by a hard sphere potential, 



V mt {r) 



oo r < a 
r > a. 



(3) 



Introducing lengths in units of the characteristic trap 
length a,ho = (fr/m tjJhn ) 1 / 2 '', r — > r/ai lQ: and energies in 
units of htuho as in |15| . the many-body Hamiltonian is: 

N 1 

H = Y1 2 ( ^ v * + r * 2) + 5 ^ (|ri - rj l} - (4) 

i i<j 



A. Diffusion Monte Carlo implementation 

Diffusion Monte Carlo is a method for finding the exact 
properties of the quantum mechanical ground-state of a 
many-body system to within arbitrary precision - see, for 
example, Ref. [55| . The starting point for this method is 
the time dependent Schroedinger equation in imaginary 
time: 



[_ v 2 + V(R) - £tMR, t) 
2m 



(5) 



The time dependent component of ^(R, t), Qi(t), is 
Qi(t) — exp[—(Ei — ET)t/h]. Et is an adjustable tar- 
get energy. In the t — > oo limit, the steady state solution 
of JSJ) is the ground state <I>o(R). 

The term diffusion Monte Carlo comes from the resem- 
blance of JSJ to the classic diffusion equation: 



DV 2 p(R,t) 



dt ' 



(6) 



This equation can be simulated by a Monte Carlo ran- 
dom walk in configuration space. Treating the [V(R) — 
Et]^(R, t) component of (0 alone results in a rate equa- 
tion of the form 



v(R)p(R, t) 



d P (R,t) 
dt ' 



(7) 



This component represents a branching process in which 
the growth or decay of a population is proportional to its 
density. In the present implementation the diffusion and 
branching processes are combined to simulate (JSJ and 
obtain the zero temperature ground state of the time 
independent Schrodinger equation. 



A simple application of (0 above results in a branching 
rate which is proportional to the potential energy U(R) — 
Et- This means that large fluctuations in the potential, 
V(R), will cause correspondingly large fluctuations in 
the population of walkers. Dramatic fluctuations in the 
number of walkers can result in large inefficiencies when 
treating realistic many-body systems. The solution to 
this problem was first presented by Kalos et al. |48| . 
In this method, a trial function is introduced to guide 
the metropolis walk to regions of higher probability and 
lower potential energy resulting in lower fluctuations in 
the population of walkers. The wave- function in JSJ is 
replaced by a product of the true ground state, ^(R, t), 
and a guiding function \I/t(R), 



*(R, t) f(R, t)* T (R)- 



(8) 



While use of a guiding function is necessary for the effi- 
cient application of the DMC method, it can introduce a 
bias into the calculation of observables which do not com- 
mute with the Hamiltonian unless corrective measures 
are taken - such as the application of "forward walking" 

mm. 

We evaluate the expectation value, (\f , |0|\I f ), of an op- 
erator O, using QMC. In integral form the expectation 
value is 

(tf = J dR **(R)C(R)*(R). (9) 

To evaluate this expression using QMC, © is recast as 

■0(R)*(R) 



= / dR |*(R 



(10) 



The result of a QMC calculation is a set of configurations 
{Ri, ...,Rm} sampled from \^>\ 2 . Using these configura- 
tions we may estimate (^ , |0|5') as 



M 



mom 



1 ^0(Ri)*(Ri) 



M 



i=i 



*(R0 



(11) 



This estimate becomes exact as M — > oo. 



B. The OBDM and natural orbitals 

A goal in this work is to describe BEC in systems with 
interactions. To do this we require a definition of the con- 
densate single particle state. Following Penrose and On- 
sager, Lowdin and others |10llll|. we take the one-body 
density matrix (OBDM) as the fundamental quantity for 
an interacting system and define the natural single par- 
ticle orbitals (NO) in terms of the OBDM. The OBDM 

is H3 



p(V,r) = <*V),*M>, 



(12) 



where ^(r) is the field operator that annihilates a single 
particle at the point r in the system. To define the NO, 
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we introduce a set of single particle states having wave 
functions <fo(r) and expand W(r) in terms of these states 
and the operators Sj which annihilate a particle from 



(13) 



directions of r' and taking the average result. In ad- 
dition, since we are dealing with identical Bosons, the 
OBDM does not depend on the particle being evaluated 
so pi(ri,r[) — jO{(rj, 7*j). This allows us to take the aver- 
age, pi(ri,r[) = l/NY,iPi(n,ri). 



Requiring that the di satisfy the usual commutation 
([ajjOj] = 5ij) and number relations ((ajaj) = NiSij), 
we have 



p(r,rO=E«^(r')^(r)^% 



(14) 



This may be taken as the defining relation of the NO, 
4>i(r). Specifically, we have from l|14|) . 



drrfr'^(r)p(r,r')^(r')=iV i 5., 



(15) 



so that the NO may be obtained by diagonalizing the 
OBDM. The eigenvectors are the NO and the eigenval- 
ues are the occupation, Ni, of the orbitals. In principle 
any orbital which satisfies Ni » 1 may be considered a 
macroscopically occupied pseudo-particle state - i.e. the 
equivalent of a Bose-Einstein condensate. A Bose sys- 
tem with more than one macroscopically occupied state 
would represent a fragmented condensate Q . In the sys- 
tems studied in this work, only a single condensate or- 
bital was found to have macroscopic occupation. The 
condensate is therefore the orbital having the highest oc- 
cupation, denoted 4>o(r), and the condensate fraction is 
n = N /N. 

The relations (|14|l and l)15|l involve the vector r and r' 
and cannot be solved directly as matrix equations. To 
obtain matrix equations, we restrict ourselves to spheri- 
cal traps and seek equations for the radial component of 
the NO as in ref. [50|. In this approach, the OBDM is 
expanded in Legendre Polynomials, Pi (f 1 • f ' ) , and eval- 
uated using the QMC ground state, ^o, a s 



ft^-O - / dn 1 dr 2 .Ar N %(r 1 ..r N )P l (r 1 -r[)Mr[--r N ). 

(16) 



C. QMC Evaluation of pi(r,r') 

In QMC we evaluate Hlt)|) in a form similar to 1|1(J|) 
giving 



Pl{n,r'x) » 



*(n,R) 



(17) 

where R = (r*2, .., tn) and e is the width of the grid ele- 
ments upon which pi(ri,r[) is being evaluated. Because 
the systems we are evaluating are spherically symmetric, 
the direction of r' is arbitrary. We may take advantage 
of this fact to reduce the statistical uncertainty in esti- 
mates of pi(r,r') by evaluating i|17|) for several different 



D. Diagonalization and error estimation 

Using the method described above, the OBDM is eval- 
uated on a grid of values of r = ie and r' = je where i 
and j are integers in the range < i,j < Q (where Q is 
a maximum cutoff). Wc may then construct the discreet 
matrix, [iep l (ie,je)je], which is readily diagonalized by 
standard matrix diagonalization methods. 

Replacing the continuous matrix pi(r,r') with the dis- 
creet matrix [iep t (ie, je)je] is a potential source of sys- 
tematic error. To avoid this problem, we evaluated each 
system with decreasing values of the grid spacing, e, such 
that e 9 +i = e q /2. The largest value of e for which no sig- 
nificant change in the calculated orbitals and occupation 
numbers occurred between e q and e g +i was then used to 
determine the condensate properties for that system. 

A second potential source of error arises in treating 
pi(r,r') (which is an infinite matrix) as a finite matrix. 
Since the trapped systems are spatially finite, the prob- 
ability of finding a particle beyond the average radius, 
R, of the cloud goes to zero very quickly. For the same 
reason, pi(r,r') ~ when either r > R or r' > R. It 
is therefore, safe to treat pi{r,r') as a finite matrix. As 
a brute force test of this assertion, we evaluated several 
systems with increasingly large cutoff values. We found 
no significant change in condensate properties calculated 
from an OBDM where r,r' < R and r, r' < 2R. 

The statistical error associated with a given orbital 
and its occupation are obtained as follows. When the 
initial OBDM, p°, is calculated the variance associated 
with each matrix element in p° is obtained. The origi- 
nal p° is assumed to represent a randomly sampled event 
from a gaussian error distribution surrounding the true 
OBDM. Based on this assertion, a set of M new OBDM's, 
{p 1 ../? }, are then generated by allowing each matrix el- 
ement to randomly vary according to its statistical error. 
Each of the new OBDM, p q , are diagonalized to obtain 
their corresponding eigenvalues, nf and eigenvectors <f>1. 
An average occupation, rij = l/M^ii', and orbital, 

4>i = fill are then obtained. The variance of these 
averages is then used as an estimate of the statistical 
error of the orbitals and occupation numbers of p . 



III. RESULTS 
A. DMC Energy 

Figure [21 shows the energy per particle calculated by 
diffusion Monte Carlo, Edmc> by variational Monte 
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FIG. 2: Diffusion Monte Carlo, Edmc, variational Monte 
Carlo, Evmco, and Gross-Pitaevskii, Ecp, energies for 
trapped hard sphere Bosons as a function of maximum den- 
sity, na 3 , in the trap. Density is varied by changing scattering 
length a, 4.3 x 10 -3 < ajaho < 0.14 where auo is the trap 
length. At higher densities Edmc clearly lies above Ecp, 
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Carlo (using the simple trial function of |50|1. Evmco, 
and using the Gross-Pitaevskii equation, E G p, of trapped 
hard core Bosons as a function of maximum density, na 3 , 
in the trap. In the dilute regime, na 3 < 10~ 4 , Edmc, 
Evmco, an d Ecp are nearly indistinguishable. The dif- 
ference in energy at na 3 — 5 x lO -5 is, for example 
10~ 3 hujho which is within the error bars of the QMC cal- 
culations. At higher densities, the DMC energy lies above 
the GP result by 3% at na 3 = lCF 3 . Evmco agrees well 
with the DMC results with a difference of only 0.3% at 
no 3 = 1CT 3 . 

Figure |31 shows the percent difference between Edmc 
and E G p, SE/E gp = (Edmc - E GP )/E GP for N = 128 
hard sphere Bosons in a spherically symmetric harmonic 
trap at higher densities, na 3 . Here, and throughout this 
paper, GP energies are calculated using a self interaction 
term proportional to (N — l)o/a^ . GP results using 
Na/aho significantly overestimate the energy for small 
N. The difference between DMC and GP energies is 
well described by 8E/Ecp oc (na 3 ) 2 / 3 . This dependence 
holds even up to trap densities of na 3 « 0.32, well above 
the density of liquid helium (na 3 k, 0.21). At this density, 
Eqp and Edmc differ by as much as 80%. 

In Fig. the dependence of 5E/Ecp on the scattering 
length, a, for N — 128 Bosons in a spherically symmetric 
harmonic trap is shown. The figure shows good agree- 
ment with SE cx (a/a^o) 8 / 5 . This is precisely the power 
law relation predicted by the first-order correction to the 
Gross-Pitaevskii energy which takes into account parti- 
cles above the condensate, denoted the modified Gross- 



FIG. 3: Percent difference between diffusion Monte Carlo, 
Edmc, and Gross-Pitaevskii, Eqp, energies for hard core 
Bosons in a spherically symmetric harmonic trap as a func- 
tion of maximum density na 3 in the trap. The % differ- 
ence between DMC and GP energies is well described by 
SE/E G p oc (na 3 ) 2/3 (dashed line). 
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FIG. 

tio of the scattering length, a, to the trap length, aho = 
(h/mcuho) 1 / 2 , for N = 128 Bosons in a spherically symmetric 
harmonic trap. The dashed line shows SE/Ecp oc (a/aho) 8 ^ 5 ■ 



Pitaevskii equation (MGP) energy |l5j . 

Figure shows the dependence of 8E = Edmc — Eqp 
on the number of particles, TV, in the trap. In this plot, 
the ratio of the scattering length to the characteristic 
length of the trap is a/ aho = 8 x apb/aho — 0.03464. 
The resulting range of densities at the center of the trap 
lie between na 3 w 8 x 10" 5 for N = 16 and na 3 ~ 6 x 
10 -4 for N — 1024. The DMC energy is approximately 
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FIG. 5: Dependence of SE — Edmc — Ecp on the num- 
ber of particles, N, in units of huiho where the ratio of the 
scattering length to the characteristic length of the trap is 
a/a ho = 0.03464. The dashed line is SE oc iV 3/5 . 
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FIG. 7: Difference between DMC and VMC0 energies 
(Evmco — Edmc) compared with the variance of the VMC0 
calculation, u(Evmco), as a function of the ratio of the hard 
sphere diameter to the trap length, a/aho- 



2% higher than the GP energy when N = 1024. The 
dashed line is a least squares fit of SE to a function of 
the form q(N) = q + qi N 3 / 5 . The relation SE(N) oc 
iV 3 / 5 is again consistent with the result obtained from 
the modified Gross-Pitaevskii equation. 



nir lQ- 5 1Q- 4 
0.15 n 1 — 



10" 



0.1 



3 
< 

g 0.05 



Fixed N 
Fixed a i 
MGP ■ 



0.05 0.1 0.15 0.2 0.25 0.3 

N 3 / r '(a/a ho f/ r ' 



0.35 



contributions of the noncondensate. If this correction is 
relevant across the entire range of systems considered, 
combining the results for the dependence of 5E on N 
and a/aho as presented in figures and [S] should provide 
a single coefficient, £, such that 



SE = ^N 3 ^(a/a ho ) 8 / 5 . 



(18) 



MGP predicts £ = 5(15) 3 / 5 /(64\/2) « 0.28. In Fig. M 
the fixed a and fixed N results are shown together along 
with the MGP prediction for the first order contribution 
to the ground state energy of atoms depleted from the 
condensate. The figure demonstrates that for systems 
with na 3 < 5 x 10~ 4 , MGP provides a good description 
of the DMC corrections to the GP energy. At higher 
densities, while the fixed a and fixed N results are sepa- 
rately well described by (a/aho) 8 ^ 5 and iV 3 / 5 respectively, 
they do not share a common coefficient £. This suggests 
that at higher densities corrections to the condensate en- 
ergy have a more complicated dependance on TV and a 
than JTBJ>- 



B. Range of validity of VMCO results 



FIG. 6: SE = {Edmc - E GP ) as a function of N 3/5 {a/a ho ) s/5 
for fixed number of particles, N = 128, (filled circles) and 
fixed scattering length, a/aho — 8 x am, (open circles) along 
with MGP prediction (heavy dashed line) . Values of the max- 
imum trap density, na 3 , for the fixed N case are shown on the 
top axis. 

The MGP expression for the ground state energy pro- 
vides the first correction to the GP energy arising from 



To investigate the range of validity of the VMCO 
trial function we evaluated the variance of the Hamil- 
tonian. If the trial function is an exact representation of 
an eigenstate of the Hamiltonian the variance is zero. 
Figure [7| provides a comparison of the difference be- 
tween DMC and VMCO results for the energy per boson, 
(EyMco — Edmc) t and the variance of the energy per 
boson, rj(EvMCo), as a function of the ratio of the hard 
sphere diameter to the trap length, a/aho- Results are 
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FIG. 8: QMC values of the width, R, of the ground 
state density of hard-core Bosons in a harmonic trap verses 
(Na/aho) 1 ^ 5 where N is the number of particles and a is the 
hard core diameter. Diamonds show the dependence when TV 
is fixed (N = 128) and a is varied, 4.33 x 10" 3 < a/a ho < 
1.11. Circles show the dependence when a is fixed (a/a^o = 
8 x a ab /a ho « 0.035) and N is varied, 32 < N < 1024. The 
short-dashed and long-dashed lines are linear and spline fits 
to the fixed a and fixed N data respectively. 



for N = 128 hard core Bosons in a spherically symmet- 
ric trap. Up to a value of a/a^o ~ 0.3, the DMC and 
VMC0 energies agree to within the variance of Evmco- 
The maximum density of the trapped Bosons for this 
"critical" value of a/a ho = 0.3 is na 3 w 3 x 10~ 2 . This 
indicates that for systems with na 3 < 10~ 2 the VMC0 
trial function not only provides a valid upper bound on 
the energy but a valid lower bound as well. 



C. Spatial distribution of trapped Bosons 

The spatial distribution of trapped Bosons is a prop- 
erty which is accessible to experimental observation. The 
first observations of BEC used the difference between a 
classical Boltzmann distribution and a condensate distri- 
bution as evidence for the existence of BEC [l2l IT3L Hif . 
Spatial resolution in most observations to date is, how- 
ever, not very high (typically only O(10 _1 ) times the 
size of the condensate itself) E3 ■ In this section, 
we compare the present QMC results for the density 
of the many-body ground state, n(r), with predictions 
of mean-field theory for the spatial distribution of the 
condensate, no(r). While this comparison is not always 
strictly correct, since depletion of the condensate means 
n(r) =/= no(r), what is actually observed in experiments is 
the "total" density which includes condensate and non- 
condensate atoms alike. The condensate distribution and 
"total" density have been treated as identical in the anal- 



ysis of experimental results |4fJ . For this reason, we will 
compare n(r) and mean field results for no(r) as if they 
are indeed measurements of the same physical quantity. 



1. The "width" of a trapped cloud of Bosons 

The radius of the condensate as predicted by the Gross- 
Pitaevskii equation in the Thomas- Fermi limit (Na >> 
1, a/a ho « 1) is 15] 



Rtf = a ho {lhN— y/ 5 . 



(19) 



We have defined the radius of the ground state, Rqmc 
by setting a cut-off value of the QMC number density, 
n(r), so that ti(Rqmc) — 10~ 5 . Figure 00 shows 
QMC results for the dependence of the width, R/aho, 
of the ground state density of N hard-core Bosons on 
the product (Na/aho) 1 ^ ■ In the figure, diamonds show 
the dependence when the number of particles is fixed, 
N = 128, and the hard core diameter, a, is varied, 
4.33 x 10~ 3 < a/a ho < 1.11. The dashed line is a 
spline fit to the fixed TV = 128 data to guide the eye. 
Circles show the dependence when the hard-core diam- 
eter is fixed, a/aho — 0.035, and N is varied, 32 < 
N < 1024. The dashed line is a linear least squares 
fit to the fixed-a data with slope w 0.52. In the region 
1 < (Na/a ho ) 1/b < 1.75, both fixed-a and fixed-TV 
results have a linear dependence on (Na/aho) 1 ^ 5 with 
the same slope. In the region where the dependence on 
(Na/aho) 1 / 5 holds, the maximum density of the trapped 



Bosons ranges from 10 



-6 < 



< 5 x 10" 



For 



small values of (Na/a,ho), the (Na/aho) 1 ^ 5 dependence 
is not expected to hold since even a single particle non- 
interacting system has a finite width. The width of the 
many-body ground state is no longer linearly dependent 
on (a/aho) 1 ^ for values of (Na/aho) 1/b > 1-75. In this 
regime, the maximum trap density is na 3 ^ 5 x 10~ 3 
and a/aho ^ 0.1. The present DMC results indicate 
that for a/aho 0.1, the width of the many-body 

ground state depends on a/aho as i? oc (a/a^ Q ) 2 / 3 rather 
than (a/aho) 1 ^ 5 ■ Linear dependence on TV 1 / 5 continues 
to hold up to the highest number of particles considered 
(N = 1024). 



2. The Total Density Profile 

Figure shows the DMC density profiles for 128 hard 
sphere trapped Bosons for four values of the maximum 
trap density na 3 . Frame (a) shows the radial density 
profile for a/aho = 8 x aab/aho — 0.03464. The max- 
imum density of the trapped Bosons in this system oc- 
curs at the center of the trap with na 3 w 2.25 x 10~ 4 . 
This corresponds to a typical density observed in ex- 
periments in metastable He* [44|> |45| . Frame (b) shows 
results for a/a^ Q = 32 x aR /aho — 0.13856 and max- 
imum density na 3 ~ 6 x 10 -3 , which is comparable 



(a) 



~i r 



(b) 




FIG. 9: DMC density profiles for hard sphere trapped 
Bosons for four values of the maximum trap density na? . 
All plots are for TV = 128 and values of a/aho ~ 
0.03464,0.13856,0.27712,1.1084 for frames (a),(b),(c), and 
(d) respectively. 



with densities found in 85 Rb experiments. In frame (c), 
a/a/io = 64 x aRb/a-ho — 0.27712. Here, local correla- 
tions in the density distribution near the center of the 
trap are readily apparent. Finally, frame (d) shows the 
density profile for a/aho — 256 x a^/aho — 1.1084. In 
this system, the hard spheres appear to have solidified in 
the center of the trap. The n(r)o 3 for this system is only 
qualitatively correct as mixed estimator bias caused by 
the guiding function used is a factor here. 



3. Comparison of DMC n(r) and Thomas - Fermi riTF(r) 

In the so called "Thomas-Fermi" (TF) approximate 
form of the Gross-Pitaevskii equation, the interaction 
term g oc Na in the GP equation is assumed to dom- 
inate the "kinetic energy" or gradient term resulting in 
an analytically solvable form of the GP equation. The TF 
approximation is expected to be valid when Na is large, 
Na >> 1, the interaction density is low, na 3 << 1, and 
the ratio of the scattering length to the characteristic 
harmonic trap length is small, a/aho << 1< 

Figure ITU1 shows a comparison of the total density dis- 
tribution calculated using DMC, «.DMc( r )j to the density 
predicted by the Thomas-Fermi approximation, 



n TF 



(r) = ((15TVa) 2/5 - x 2 )/8nNa. 



(20) 




The top frame, (a), shows the density profile for N = 
1024 Bosons with a/a ho = 8x a Rb = 0.03464. Here, the 



FIG. 10: Comparison of total density distribution calcu- 
lated using diffusion Monte Carlo, riDMc('"), for hard sphere 
Bosons in a harmonic trap to the density predicted by the 
Thomas- Fermi approximation 12011 . Top frame (a) is for 
TV = 1024 Bosons with ratio of scattering length to trap 
length of a/aho = 8 x ajib = 0.03464. Density is expressed in 
terms of n(r)a 3 x 10 4 . Frame (b) is for TV = 128 Bosons with 
a/fflfco = 64 x ajjb = 0.27712. Density is expressed in terms of 
n(r)a 3 x 10 2 . 



TF and DMC results agree quite well. However, the TF 
result slightly overestimates the density near the center of 
the trap and fails to reproduce the low density tail which 
occurs near the edge of the trapped cloud. Frame (b) 
shows n(r) for TV = 128 and a/a ho — 64 x am — 0.27712. 
Note that the product, Na/ai lo , (the only variable re- 
sponsible for determining the shape of the TF and GP 
density profiles) is the same in both frames. Clearly, in 
the bottom frame, the TF approximation dramatically 
overestimates the density at the center of the trap and 
underestimates the width of the condensate. 
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FIG. 11: Condensate fraction, no, as a function of the density, 
no 3 , for N = 128 trapped hard sphere Bosons. Here n is the 
number density at the center of the trap and a is the scattering 
length. Circles are from the mean-field Bogoliubov (MFB) 
expression for no in a uniform dilute Bose gas integrated over 
the TF density. The up and down-facing triangles are the 
DMC and VMC results respectfully. Dashed lines are spline 
fits to guide the eye. 



D. Condensate fraction 

FigurelrDshows the condensate fraction, no, as a func- 
tion of the density na 3 for N = 128 trapped hard sphere 
Bosons. The n is the maximum number density which is 
at the center of the trap in the density range shown. The 
density was varied by changing the value of a/aho- The 
corresponding values of a/aho are shown on the top axis. 
Circles are the mean-field Bogoliubov (MFB) result for a 
uniform dilute Bose gas integrated over the GP density 
in the Thomas-Fermi limit (5a | obtained by solving 



n = 1 - O.S798(N a/a ho ) 6/5 /N. 



(21) 



The up and down-facing triangles are the DMC and 
VMCO results, respectively, obtained from diagonalizing 
the OBDM. For na 3 < 10 -4 , all three values of no agree 
to within 1%. At higher densities, the MFB result consis- 
tently overestimates the condensate fraction. MFB over- 
estimates the condensate fraction because it ignores lo- 
cal pair correlations which act to deplete the condensate. 
In contrast the DMC value of the condensate is consis- 
tently lower than either the VMCO or MFB estimates. 
We believe that the DMC result for no is lower than ei- 
ther VMC and MFB because it is able to treat local pair 
correlations more accurately. Pair correlations allow the 
total energy to decrease at the expense of long range or- 
der. Since DMC is able to sample the exact ground state, 
the mixed estimate for no obtained from DMC is more 
accurate than VMC or MFB. 
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FIG. 12: Condensate fraction, no, over a wide density range. 
The legend is the same as Fig. 1111 



TABLE I: Condensate fraction as obtained from mean-field- 
Bogoliubov (MFB) , VMC, and DMC methods. 



na 3 




a/aho 


MFB 


VMC 


DMC 


1.3 x 10" 


-6 


1 


0.998 


0.999(9) 


0.99(9) 


4.6 x 10" 


-5 


4 


0.992 


0.992(4) 


0.99(2) 


2.5 x 10" 


-4 


8 


0.983 


0.977(8) 


0.97(7) 


2.5 x 10" 


-3 


16 


0.959 


0.942(1) 


0.94(2) 


2.4 x 10" 


-2 


64 


0.785 


0.745(7) 


0.70(2) 


1.1 x 10" 


-1 


128 


0.506 


0.476(5) 


0.3(5) 


3.2 x 10" 


-1 


256 


N/A a 


0.160(0) 


0.1(0) 



"MFB predicts a negative condensate fraction for this system. 



Figure 1121 shows no over a wider density range. Here 
n is again the maximum number density in the trap. 
At high densities the maximum density in the trap is 
not always at the center of the trap. As in the dilute 
regime presented in Fig. ^1 MFB consistently overes- 
timates the condensate fraction for most densities. At 
na 3 s» 0.28, however, the MFB estimate of no goes to 
zero while both VMC and DMC still show a condensate 
fraction of no w 10%. The MFB estimate goes to zero be- 
cause the TF density profile used to calculate the MFB 
value of no does not have a broad low density region 
near the surface of the trapped cloud of atoms as do the 
DMC and VMC density distributions. As will be shown 
in the next section, the dilute region at the "edge" of the 
trapped cloud can support a condensate even when the 
condensate fraction in the dense center of the trap goes 
to zero. 

The density corresponding to liquid helium at SVP 
(na 3 — 0.21) is indicated on the plot. At this density, 
VMC gives a condensate fraction of no ~ 25% while 
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E. Spatially dependent depletion of the condensate 

In Fig. ^|we compare the total density distribution, 
n(r), to the condensate distribution, ng{r) = no|</>o( r )| 2 , 
for N = 128 hard sphere Bosons in a harmonic trap cal- 
culated using diffusion Monte Carlo. In the top frame, 
a/a-ho — 64 x anb/aho giving a maximum density in the 
trap of na 3 w 2.4 x 10~ 2 and total condensate fraction 
of no w 70%. In this system, the spatial distribution 
of the condensate follows the shape of the total density 
distribution except at small r. It is worth noting that 
while the total density exhibits local correlations in the 
dense region near the center of the trap, the condensate 
distribution is relatively flat in this region. In the bot- 
tom frame of the figure, a/ai lo — 256 x am/O'ho = 1-1084 
resulting in a maximum density of na 3 m 0.325 and a 
condensate fraction of no w 10%. This is the same sys- 
tem shown in frame (d) of Fig. |5J As discussed above, 
the DMC results at this density are biased by the VMC 
guiding function used. Nevertheless, we believe the re- 
sults to be qualitatively correct. Here, strong pair cor- 
relations have completely depleted the condensate in the 
center of the trap but the relatively dilute region near 
the edge of the trap is still able to support a conden- 
sate. We find that for trapped hard sphere Bosons, the 
local condensate fraction, n (r)/n(r), rises in the dilute 
region near the surface and remains close to one all the 
way to the surface of the cloud. This may be contrasted 
with predictions for no{f)/n(r) for self-bound superfluid 
He at a free surface in which surface correlations sig- 
nificantly _dep_lete the condensate at the liquid-vacuum 
interface 



IV. DISCUSSION 
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FIG. 13: Comparison of total density distribution, n(r), to 
condensate distribution, no(r) = \4>(r)\ 2 , for TV = 128 hard 
sphere Bosons in a harmonic trap calculated using diffusion 
Monte Carlo. Circles are the total density while triangles 
represent the condensate. Dashed lines are spline fits to guide 
the eye. In the top frame, the maximum density in the trap 
is na 3 m 2.4 x 10 -2 and the total condensate fraction is no » 
70%. In the bottom frame, na 3 « 2.4 x 10~ 2 and n » 10%. 



DMC estimates a condensate fraction of uq s» 18%. 
In bulk liquid 4 He, the condensate fraction is no rs 
7.25% ■8]. This difference is explained by the fact that the 
dilute region near the surface of the trapped cloud allows 
for a larger fraction of particles to occupy the condensate 
orbital than in an uniform system at 4 He densities. 



TableUsummarizes the present DMC and VMC results 
for the condensate fraction over a wide density range. 



The main objectives of this work are to explore the 
role of interactions in determining the zero temperature 
properties of the trapped Bose gas over a wide range of 
densities and to determine the limits of the mean-field 
description of the condensate properties. To this end, we 
have employed quantum Monte-Carlo (QMC) methods 
and the one-body density matrix (OBDM) formulation 
of BEC. We find the OBDM description of a many-body 
Bose system combined with QMC techniques, provides a 
coherent method for the study of the ground state prop- 
erties and Bose-Einstein condensation in traps from the 
dilute to the very dense regime. By comparing our QMC 
results with mean-field theory we determine key limits of 
the mean field description. 



A. The ground state energy 



We find that in the dilute limit, na 3 < 10~ 4 , where 
the condensate depletion is small, no > 99%, the GP de- 
scription of the condensate provides a good description 
of the full many-body ground state. Once the density has 



reached, na 3 w 10 -3 , approximately 6% of the atoms lie 
outside of the condensate and the condensate energy ob- 
tained from GP theory lies 3% below the QMC energy. 
For na 3 > 1CP 3 , the GP energy does not describe the 
energy of the Bosons in the trap accurately. The present 
QMC corrections to the GP energy, SE = Edmg — Egp, 
are proportional to N 3 / 5 when N is allowed to vary with 
fixed a and are proportional to (a/ciho) 8 ^ 5 when a is al- 
lowed to vary with fixed N. This dependence on N and 
a holds for all densities studied (10 -6 < na 3 < 0.5) 
and is consistent with the expected corrections to the 
GP energy arising from the depletion of the conden- 
sate (JTSJ. Thus, the GP description of the condensate 
energy appears to be valid even in the highly interacting 
regime. However, the dependence of SE on the prod- 
uct N 3 / 5 (a/a ho ) 8/5 , as predicted by MGP (JTHJ, holds 
up to densities na 3 ~ 5 x 10~ 4 only. As interaction is 
increased the effects of the non-condensate play an in- 
creasingly significant roll in determining the properties 
of the total ground state and a more complicated func- 
tional dependence of 5E(N, a) than the simple product 
iV 3 / 5 (a/a; lo ) 8 / 5 is required at higher densities. 



B. Deviations from the mean-field description 

Figure IT4l contains striped bands indicating regions in 
a/o-ho where QMC results diverge from mean-field / Bo- 
goliubov predictions. Since the degree of depletion of the 
condensate arising from inter-Boson interaction plays a 
significant role in determining beyond-mean-ficld effects, 
QMC values for the number of atoms outside the con- 
densate, N, for a system with a total of N — 128 Bosons 
are shown along with the regions. The first sign of di- 
vergence (a) occurs at a density of na 3 w 3 x 10~ 4 and 
a value of a/aho ~ 8 x anb/o-ho ~ 0.035. At this density, 
QMC and MFB (|21 I) results for the condensate fraction, 
?io = No/N, begin to diverge (see Fig. EJ). Below this 
value of a/aho, QMC and MFB values of uq agree to 
within 1%. At na 3 ps 10~ 3 , MFB underestimates the de- 
pletion of the condensate by 30%. At higher densities, 
na 3 ps 10" 1 , MFB predicts a condensate fraction 40% 
higher than QMC. 

The second point of interest in Fig. [21 marked (b) 
occurs in the region of na 3 ps 2.5x 10~ 3 and a/aho ~ 0.12. 
Near this value of a/aho, QMC results for the size of the 
many-body ground state and mean-field results for the 
size of the condensate (|19|l begin to differ. For values of 
a/aho ^ 0.12, the width of the many-body ground state 
is proportional to (Na/aho) 1 ^ 5 as predicted by mean- field 
theory. At higher values, (a/aho > 0.12), we find that the 
size of the condensate is better described by a scaling of 
(a/aho) 2 / 3 ■ The scaling is shown in Fig. [S|and discussed 
in Section fill C II Thus, for systems with na 3 > 10~ 3 , 
GP theory in the TF limit under-estimates the growth of 
the size of the ground state with a/a} lo significantly. In 
the extreme range of very large scattering length or very 
tight trapping potential where a/aho = 1, GP predicts 
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FIG. 14: QMC determination of the number of Bosons out- 
side the condensate, N, for TV = 128 hard sphere Bosons 
in a spherically symmetric harmonic trap vs. a/aho- The 
corresponding densities, na 3 , are shown across the top axis. 
Striped bands indicate regions in a/aho where the QMC re- 
sults for BEC properties diverge from mean-field / Bogoliubov 
predictions, (a) QMC and Bogoliubov results for no begin to 
diverge, (b) QMC and mean-field results for the size of the 
condensate diverge, (c) Local correlations in the density pro- 
file of the many-body ground state begin to appear, (d) Con- 
densate begins to shift from center of trap, (e) Condensate 
exists only in dilute region near surface of trapped cloud. 



a condensate distribution 20% smaller than the width of 
the ground state obtained from QMC. 

The band (c) in Fig. [21 indicates the region in which 
local correlations in the density profile of the many- 
body ground state begin to appear. These local corre- 
lations signal a clear departure from mean-field proper- 
ties. This effect occurs for systems with trap densities of 
na 3 > 2.5 x 10~ 2 and a/a ho w 64 x a Rb /aho ~ 0.28. 
At this level of interaction the condensate fraction as ob- 
tained from DMC is uq ps 70%. We find that at this den- 
sity the condensate density is smoothly varying through- 
out the trap with little or no local density fluctuations 
(See top frame of Fig. I13|) . Evidence that the condensate 
distribution does not explicitly follow the total density 
distribution is another demonstration that a local den- 
sity approximation description of the condensate breaks 
down at this density. 

The band marked (d) in Fig. 1141 approximates the re- 
gion in which the condensate begins to shift from the 
center of the trap to the surface. Here, na 3 ps 0.2 and 
a/aho ~ 0.8. The condensate fraction is n ~ 20%. At 
this level of interaction and beyond, mean-field approxi- 
mations and the Bogoliubov approximation both fail to 
appropriately describe the properties of a trapped BEC. 
We speculate that at this density, the increased deple- 
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tion in the center of the trap could effectively pin vortex 
states. 

The final point of interest in Fig. ll4l occurs in the region 
marked by the band (e). For systems with na 3 > 0.3 and 
fl / a /io ^ 256 X <XRb/aho w 1-1) the condensate exists only 
in dilute region near surface of trapped cloud. Strong pair 
correlations have completely depleted the condensate in 
the center of the trap but the relatively dilute region near 
the edge of the trap is still able to support a condensate. 
Figure ITTA presents DMC results which demonstrate this 
phenomena. 
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